Delphi’s COVIDcast Project:
An Ecosystem for Tracking and Forecasting the Pandemic
October 9, 2020
October 9, 2020
How many people have died from COVID-19 per day, in my state, since March 1?
library(covidcast)
deaths = covidcast_signal(data_source = "usa-facts",
signal = "deaths_7dav_incidence_num",
start_day = "2020-03-01", end_day = "2020-09-25",
geo_type = "state", geo_values = "pa")
plot(deaths, plot_type = "line",
title = "COVID-19 deaths in PA (7-day trailing average)")
What percentage of daily hospital admissions are due to COVID-19 in PA, NY, TX?
hosp = covidcast_signal(data_source = "hospital-admissions",
signal = "smoothed_adj_covid19",
start_day = "2020-03-01", end_day = "2020-09-25",
geo_type = "state", geo_values = c("pa", "ny", "tx"))
plot(hosp, plot_type = "line",
title = "% of hospital admissions due to COVID-19")
What does the current COVID-19 incident case rate look like, nationwide?
cases = covidcast_signal(data_source = "usa-facts",
signal = "confirmed_7dav_incidence_prop",
start_day = "2020-09-25", end_day = "2020-09-25")
plot(cases, title = "Daily new COVID-19 cases per 100,000 people")
What does the current COVID-19 cumulative case rate look like, nationwide?
cases = covidcast_signal(data_source = "usa-facts",
signal = "confirmed_cumulative_prop",
start_day = "2020-09-25", end_day = "2020-09-25")
plot(cases, title = "Cumulative COVID-19 cases per 100,000 people",
choro_params = list(legend_n = 6))
Where is the current COVID-19 cumulative case rate greater than 2%?
plot(cases, choro_col = c("#D3D3D3", "#FFC0CB"),
title = "Cumulative COVID-19 cases per 100,000 people",
choro_params = list(breaks = c(0, 2000), legend_width = 3))
How do some cities compare in terms of doctor’s visits due to COVID-like illness?
dv = covidcast_signal(data_source = "doctor-visits",
signal = "smoothed_adj_cli",
start_day = "2020-03-01", end_day = "2020-09-25",
geo_type = "msa",
geo_values = name_to_cbsa(c("Pittsburgh", "New York",
"San Antonio", "Miami")))
plot(dv, plot_type = "line",
title = "% of doctor's visits due to COVID-like illness")
How do my county and my friend’s county compare in terms of COVID symptoms?
sympt = covidcast_signal(data_source = "fb-survey",
signal = "smoothed_hh_cmnty_cli",
start_day = "2020-04-15", end_day = "2020-09-25",
geo_values = c(name_to_fips("Allegheny"),
name_to_fips("Fulton", state = "GA")))
plot(sympt, plot_type = "line", range = range(sympt$value),
title = "% of people who know somebody with COVID symptoms")
The last two weeks of August in CA …
# Let's get the data that was available as of 09/22, for the end of August in CA
dv = covidcast_signal(data_source = "doctor-visits",
signal = "smoothed_adj_cli",
start_day = "2020-08-15", end_day = "2020-08-31",
geo_type = "state", geo_values = "ca",
as_of = "2020-09-21")
# Plot the time series curve
xlim = c(as.Date("2020-08-15"), as.Date("2020-09-21"))
ylim = c(3.83, 5.92)
ggplot(dv, aes(x = time_value, y = value)) +
geom_line() +
coord_cartesian(xlim = xlim, ylim = ylim) +
geom_vline(aes(xintercept = as.Date("2020-09-21")), lty = 2) +
labs(color = "as of", x = "Date", y = "% doctor's visits due to CLI in CA") +
theme_bw() + theme(legend.pos = "bottom")
The last two weeks of August in CA …
# Now loop over a bunhch of "as of" dates, fetch data from the API for each one
as_ofs = seq(as.Date("2020-09-01"), as.Date("2020-09-21"), by = "week")
dv_as_of = map_dfr(as_ofs, function(as_of) {
covidcast_signal(data_source = "doctor-visits", signal = "smoothed_adj_cli",
start_day = "2020-08-15", end_day = "2020-08-31",
geo_type = "state", geo_values = "ca", as_of = as_of)
})
# Plot the time series curve "as of" September 1
dv_as_of %>%
filter(issue == as.Date("2020-09-01")) %>%
ggplot(aes(x = time_value, y = value)) +
geom_line(aes(color = factor(issue))) +
geom_vline(aes(color = factor(issue), xintercept = issue), lty = 2) +
coord_cartesian(xlim = xlim, ylim = ylim) +
labs(color = "as of", x = "Date", y = "% doctor's visits due to CLI in CA") +
geom_line(data = dv, aes(x = time_value, y = value)) +
geom_vline(aes(xintercept = as.Date("2020-09-21")), lty = 2) +
theme_bw() + theme(legend.pos = "none")
The last two weeks of August in CA …
dv_as_of %>%
ggplot(aes(x = time_value, y = value)) +
geom_line(aes(color = factor(issue))) +
geom_vline(aes(color = factor(issue), xintercept = issue), lty = 2) +
coord_cartesian(xlim = xlim, ylim = ylim) +
labs(color = "as of", x = "Date", y = "% doctor's visits due to CLI in CA") +
geom_line(data = dv, aes(x = time_value, y = value)) +
geom_vline(aes(xintercept = as.Date("2020-09-21")), lty = 2) +
theme_bw() + theme(legend.pos = "none")
Why ask a proxy question (have people report on others)? Here’s Spearman correlations to COVID-19 case rates sliced by time:
# Fetch Facebook % CLI signal, % CLI-in-community signal and confirmed case
# incidence proportions
start_day = "2020-04-15"
end_day = "2020-09-05"
sympt1 = covidcast_signal("fb-survey", "smoothed_cli",
start_day, end_day)
sympt2 = covidcast_signal("fb-survey", "smoothed_hh_cmnty_cli",
start_day, end_day)
cases = covidcast_signal("usa-facts", "confirmed_7dav_incidence_prop",
start_day, end_day)
# Consider only counties with at least 500 cumulative cases so far
case_num = 500
geo_values = covidcast_signal("usa-facts", "confirmed_cumulative_num",
max(cases$time), max(cases$time)) %>%
filter(value >= case_num) %>% pull(geo_value)
sympt1_act = sympt1 %>% filter(geo_value %in% geo_values)
sympt2_act = sympt2 %>% filter(geo_value %in% geo_values)
cases_act = cases %>% filter(geo_value %in% geo_values)
# Compute correlations, per time, over all counties
df_cor1 = covidcast_cor(sympt1_act, cases_act, by = "time_value",
method = "spearman")
df_cor2 = covidcast_cor(sympt2_act, cases_act, by = "time_value",
method = "spearman")
# Stack rowwise into one data frame
df_cor = rbind(df_cor1, df_cor2)
df_cor$signal = c(rep("% CLI", nrow(df_cor1)),
rep("% CLI-in-community", nrow(df_cor2)))
# Then plot correlations over time
ggplot_colors = c("#FC4E07", "#00AFBB", "#E7B800")
ggplot(df_cor, aes(x = time_value, y = value)) +
geom_line(aes(color = signal)) +
scale_color_manual(values = ggplot_colors[c(3,1)]) +
labs(title = "Correlation between survey signals and case rates (by time)",
subtitle = sprintf("Over all counties with at least %i cumulative cases",
case_num), x = "Date", y = "Correlation") +
theme_bw() + theme(legend.pos = "bottom", legend.title = element_blank())
Now here’s Spearman correlations to COVID-19 case rates sliced by county:
# Compute correlations, per time, over all counties
df_cor1 = covidcast_cor(sympt1_act, cases_act, by = "geo_value",
method = "spearman")
df_cor2 = covidcast_cor(sympt2_act, cases_act, by = "geo_value",
method = "spearman")
# Stack rowwise into one data frame
df_cor = rbind(df_cor1, df_cor2)
df_cor$signal = c(rep("% CLI", nrow(df_cor1)),
rep("% CLI-in-community", nrow(df_cor2)))
# Then plot correlations as densities
ggplot(df_cor, aes(value)) + geom_density(aes(color = signal, fill = signal),
alpha = 0.4) +
scale_color_manual(values = ggplot_colors[c(3,1)]) +
scale_fill_manual(values = ggplot_colors[c(3,1)]) +
labs(title = "Correlation between survey signals and case rates (by county)",
subtitle = sprintf("Over all counties with at least %i cumulative cases",
case_num), x = "Date", y = "Correlation") +
theme_bw() + theme(legend.pos = "bottom", legend.title = element_blank())
Let’s take a look at case counts in Miami-Dade, from June 1 to July 15, and compare it to the % CLI-in-community indicator based on our survey:
# Fetch Facebook % CLI-in-community signal and confirmed case incidence numbers
# from June 1 to July 15
start_day = "2020-06-01"
end_day = "2020-07-15"
sympt = covidcast_signal("fb-survey", "smoothed_hh_cmnty_cli",
start_day, end_day)
cases = covidcast_signal("usa-facts", "confirmed_7dav_incidence_num",
start_day, end_day)
# Function to transform from one range to another
trans = function(x, from_range, to_range) {
(x - from_range[1]) / (from_range[2] - from_range[1]) *
(to_range[2] - to_range[1]) + to_range[1]
}
# Function to produce a plot comparing the signals for one county
plot_one = function(geo_value, title = NULL, xlab = NULL,
ylab1 = NULL, ylab2 = NULL, legend = TRUE) {
# Filter down the signal data frames
given_geo_value = geo_value
sympt_one = sympt %>% filter(geo_value == given_geo_value)
cases_one = cases %>% filter(geo_value == given_geo_value)
# Compute ranges of the two signals
range1 = cases_one %>% select("value") %>% range
range2 = sympt_one %>% select("value") %>% range
# Convenience functions for our two signal ranges
trans12 = function(x) trans(x, range1, range2)
trans21 = function(x) trans(x, range2, range1)
# Find state name, find abbreviation, then set title
state_name = fips_to_name(paste0(substr(geo_value, 1, 2), "000"))
state_abbr = name_to_abbr(state_name)
title = paste0(fips_to_name(geo_value), ", ", state_abbr)
# Transform the combined signal to the incidence range, then stack
# these rowwise into one data frame
df = select(rbind(sympt_one %>% mutate_at("value", trans21),
cases_one), c("time_value", "value"))
df$signal = c(rep("% CLI-in-community", nrow(sympt_one)),
rep("New COVID-19 cases", nrow(cases_one)))
# Finally, plot both signals
pos = ifelse(legend, "bottom", "none")
return(ggplot(df, aes(x = time_value, y = value)) +
geom_line(aes(color = signal)) +
scale_color_manual(values = ggplot_colors[1:2]) +
scale_y_continuous(name = ylab1, limits = range1,
sec.axis = sec_axis(trans = trans12,
name = ylab2)) +
labs(title = title, x = xlab) + theme_bw() +
theme(legend.pos = pos, legend.title = element_blank()))
}
# Produce a plot for Miami-Dade, and add vertical lines
plot_one(name_to_fips("Miami-Dade"), xlab = "Date",
ylab1 = "Daily new confirmed COVID-19 cases",
ylab2 = "% of people who know someone with CLI") +
geom_vline(xintercept = as.numeric(as.Date("2020-06-19")),
linetype = 2, size = 1, color = ggplot_colors[1]) +
geom_vline(xintercept = as.numeric(as.Date("2020-06-25")),
linetype = 2, size = 1, color = ggplot_colors[2])
Ok, that was just one county… let’s look at the top 20 in terms of the rise in case counts:
num = 20
geo_values = cases %>% group_by(geo_value) %>%
summarize(diff = last(value) - first(value)) %>%
arrange(desc(diff)) %>% head(num) %>% pull(geo_value)
p_list = vector("list", num)
for (i in 1:num) {
p_list[[i]] = plot_one(geo_values[i], legend = FALSE)
}
do.call(grid.arrange, c(p_list, nrow = 5, ncol = 4))
Results from forecasts made over early May to late August (for details, read this blog post):
# This RData file was downloaded from https://github.com/cmu-delphi/delphi-blog/tree/main/content/post/forecast-demo;
# the code for generating this RData file is also there
load("demo-extended.rda")
# Compute and plot median errors as function of number of days ahead
err_by_lead = res %>%
select(-c(err3, err4)) %>%
drop_na() %>% # Restrict to common time
mutate(err1 = err1 / err0, err2 = err2 / err0) %>% # Compute relative error
# to strawman model
ungroup() %>%
select(-err0) %>%
pivot_longer(names_to = "model", values_to = "err",
cols = -c(geo_value, time_value, lead)) %>%
mutate(model = factor(model, labels = c("Cases", "Cases + Facebook"))) %>%
group_by(model, lead) %>%
summarize(err = median(err)) %>%
ungroup()
ggplot(err_by_lead, aes(x = lead, y = err)) +
geom_line(aes(color = model)) +
geom_point(aes(color = model)) +
geom_hline(yintercept = err_by_lead %>%
filter(lead %in% 7, model == "Cases") %>% pull(err),
linetype = 2, color = "gray") +
scale_color_manual(values = c("black", ggplot_colors[1])) +
labs(title = "Forecasting errors by number of days ahead",
subtitle = sprintf("Over all counties with at least %i cumulative cases",
case_num),
x = "Number of days ahead", y = "Median scaled error") +
theme_bw() + theme(legend.pos = "bottom", legend.title = element_blank())